!/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!*  Copyleft (!c) 2004-2005 Div. Nucl. Med., Cyric, Tohoku Univ.    *
!*  Copyleft (!c) 2005-2011 Elseif Laboratory.                      *
!*  elseifkk@users.sf.net                                           *
!*                                                                  *
!*  All Rights Reversed.                                            *
!*                                                                  *
!*                                                                  *
!*  LICENSE AGREEMENT:                                              *
!*                                                                  *
!*  You agree that we don't agree to anything.                      *
!*                                                                  *
!*                                                                  *
!*  DISCLAIMER OF WARRANTIES:                                       *
!*                                                                  *
!*  This software is provided to you "AS IS," and the authors and   *
!*  contributors disclaim any warranty or liability obligations to  *
!*  you of any kind, implied or express, including but not limited  *
!*  to the warranties of merchantability, fitness for a particular  *
!*  purpose, and non-infringement.                                  *
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
!
! math function utility.
!
! TTD: 
!   none.
! note: 
!   1) some procs were written in asm as kutyA.asm in order to 
!      reduce computation time.
!
module kuty
  implicit none
  interface sort
     module procedure sort_r8
     module procedure sort_i4
     module procedure sort_i1
     module procedure sort_i2
     module procedure sort_r4
  end interface
  interface persort
     module procedure persort_r8
     module procedure persort_i4
  end interface
  private pivot_i4
  private partition_i4
  private quickSort_i4
  private pivot_r8
  private partition_r8
  private quickSort_r8
  private partitionP_r8
  private quickSortP_r8
  integer*4,pointer,private::permutation(:)
  integer*4,pointer,private::sorted(:)
contains
  !------------------------------------------------------!
  ! sorting for integer
  !------------------------------------------------------!
  integer*4 function pivot_i4(a,i,j)
    ! internal use only
    integer*4,intent(in)::a(:)
    integer*4,intent(in)::i,j
    integer*4 k
    k=i+1
    do while(k<=j.and.a(i)==a(k))
       k=k+1
    end do
    if(k>j) then
       pivot_i4=-1
    elseif(a(i)>=a(k)) then
       pivot_i4=i
    else
       pivot_i4=k
    endif
  end function pivot_i4
  !------------------------------------------------------!
  integer*4 function partition_i4(a,i,j,x)
    integer*4,intent(inout)::a(:)
    integer*4,intent(in)::i,j
    integer*4,intent(in)::x
    integer*4 l,r
    integer*4 t
    l=i
    r=j
    do while(l<=r)
       do while(l<=j.and.a(l)<x)
          l=l+1
       enddo
       do while(r>=i.and.a(r)>=x) 
          r=r-1
       enddo
       if(l>r) exit
       t=a(l)
       a(l)=a(r)
       a(r)=t
       l=l+1
       r=r-1
    end do
    partition_i4=l
  end function partition_i4
  !------------------------------------------------------!
  recursive subroutine quickSort_i4(a,i,j)
    integer*4,intent(inout),target::a(:)
    integer*4,intent(in)::i,j
    interface
       function pivot_i4_asm(a,i,j)
         integer*4,intent(inout)::a
         integer*4,intent(in)::i
         integer*4,intent(in)::j
       end function pivot_i4_asm
    end interface
    integer*4 k
    k=pivot_i4_asm(a(1),i,j)
    if(k==0) return
    if(k-1>i) call quickSort_i4(a,i,k-1)
    if(k<j) call quickSort_i4(a,k,j)
    sorted=>a
  end subroutine quickSort_i4
  !------------------------------------------------------!
  integer*4 function partitionP_i4(a,i,j,x)
    integer*4,intent(inout)::a(:)
    integer*4,intent(in)::i,j
    integer*4,intent(in)::x
    integer*4 l,r,k
    integer*4 t
    l=i
    r=j
    do while(l<=r)
       do while(l<=j.and.a(l)<x)
          l=l+1
       enddo
       do while(r>=i.and.a(r)>=x) 
          r=r-1
       enddo
       if(l>r) exit
       t=a(l)
       a(l)=a(r)
       a(r)=t
       k=permutation(l)
       permutation(l)=permutation(r)
       permutation(r)=k
       l=l+1
       r=r-1
    end do
    partitionP_i4=l
  end function partitionP_i4
  !------------------------------------------------------!
  recursive subroutine quickSortP_i4(a,i,j)
    integer*4,intent(inout),target::a(:)
    integer*4,intent(in)::i,j
    integer*4 p,k
    integer*4 x
    if(i==j) return
    p=pivot_i4(a,i,j)
    if(p/=-1) then
       x=a(p)
       k=partitionP_i4(a,i,j,x)
       call quickSortP_i4(a,i,k-1)
       call quickSortP_i4(a,k,j)
    end if
  end subroutine quickSortP_i4
  !------------------------------------------------------!
  subroutine persort_i4(n,a,p)
    integer*4,intent(in)::n
    integer*4,intent(inout)::a(n)
    integer*4,intent(out),target::p(n)
    integer*4 i
    do i=1,n
       p(i)=i
    end do
    permutation=>p
    call quickSortP_i4(a,1,n)
    nullify(permutation)
  end subroutine persort_i4
  !
  subroutine sort_i4(n,a)
    integer*4,intent(in)::n
    integer*4,intent(inout)::a(n)
    call quickSort_i4(a,1,n)
  end subroutine sort_i4
  !
  subroutine sort_i2(n,a)
    integer*4,intent(in)::n
    integer*2,intent(inout)::a(n)
    integer*4 int_a(n)
    call quickSort_i4(int_a,1,n)
    a=sorted
  end subroutine sort_i2
  !
  subroutine sort_i1(n,a)
    integer*4,intent(in)::n
    integer*1,intent(inout)::a(n)
    integer*4 int_a(n)
    call quickSort_i4(int_a,1,n)
    a=sorted
  end subroutine sort_i1
  !------------------------------------------------------!
  ! sorting for real
  !------------------------------------------------------!
  integer*4 function pivot_r8(a,i,j)
    real*8,intent(in)::a(:)
    integer*4,intent(in)::i,j
    integer*4 k
    k=i+1
    do while(k<=j.and.a(i)==a(k))
       k=k+1
    end do
    if(k>j) then
       pivot_r8=-1
    elseif(a(i)>=a(k)) then
       pivot_r8=i
    else
       pivot_r8=k
    endif
  end function pivot_r8
  !------------------------------------------------------!
  integer*4 function partition_r8(a,i,j,x)
    real*8,intent(inout)::a(:)
    integer*4,intent(in)::i,j
    real*8,intent(in)::x
    integer*4 l,r
    real*8 t
    l=i
    r=j
    do while(l<=r)
       do while(l<=j.and.a(l)<x)
          l=l+1
       enddo
       do while(r>=i.and.a(r)>=x) 
          r=r-1
       enddo
       if(l>r) exit
       t=a(l)
       a(l)=a(r)
       a(r)=t
       l=l+1
       r=r-1
    end do
    partition_r8=l
  end function partition_r8
  !------------------------------------------------------!
  recursive subroutine quickSort_r8(a,i,j)
    real*8,intent(inout)::a(:)
    integer*4,intent(in)::i,j
    integer*4 k
    interface 
       function pivot_r8_asm(a,i,j)
         real*8,intent(inout)::a
         integer*4,intent(in)::i
         integer*4,intent(in)::j
       end function pivot_r8_asm
    end interface
    k=pivot_r8_asm(a(1),i,j)
    if(k==0) return
    if(k-1>i) call quickSort_r8(a,i,k-1)
    if(k<j) call quickSort_r8(a,k,j)
  end subroutine quickSort_r8
  !------------------------------------------------------!
  subroutine sort_r8(n,a)
    integer*4,intent(in)::n
    real*8,intent(inout)::a(n)
    if(n==1) return
    call quickSort_r8(a,1,n)
  end subroutine sort_r8
  !------------------------------------------------------!
  recursive subroutine quickSort_r4(a,i,j)
    real*4,intent(inout)::a(:)
    integer*4,intent(in)::i,j
    integer*4 k
    interface 
       function pivot_r4_asm(a,i,j)
         real*4,intent(inout)::a
         integer*4,intent(in)::i
         integer*4,intent(in)::j
       end function pivot_r4_asm
    end interface
    k=pivot_r4_asm(a(1),i,j)
    if(k==0) return
    if(k-1>i) call quickSort_r4(a,i,k-1)
    if(k<j) call quickSort_r4(a,k,j)
  end subroutine quickSort_r4
  !------------------------------------------------------!
  subroutine sort_r4(n,a)
    integer*4,intent(in)::n
    real*4,intent(inout)::a(n)
    if(n==1) return
    call quickSort_r4(a,1,n)
  end subroutine sort_r4
  !------------------------------------------------------!
  ! sorting with permutation
  !------------------------------------------------------!
  integer*4 function partitionP_r8(a,i,j,x)
    real*8,intent(inout)::a(:)
    integer*4,intent(in)::i,j
    real*8,intent(in)::x
    integer*4 l,r,k
    real*8 t
    l=i
    r=j
    do while(l<=r)
       do while(l<=j.and.a(l)<x)
          l=l+1
       enddo
       do while(r>=i.and.a(r)>=x) 
          r=r-1
       enddo
       if(l>r) exit
       t=a(l)
       a(l)=a(r)
       a(r)=t
       k=permutation(l)
       permutation(l)=permutation(r)
       permutation(r)=k
       l=l+1
       r=r-1
    end do
    partitionP_r8=l
  end function partitionP_r8
  !------------------------------------------------------!
  recursive subroutine quickSortP_r8(a,i,j)
    real*8,intent(inout)::a(:)
    integer*4,intent(in)::i,j
    integer*4 p,k
    real*8 x
    if(i==j) return
    p=pivot_r8(a,i,j)
    if(p/=-1) then
       x=a(p)
       k=partitionP_r8(a,i,j,x)
       call quickSortP_r8(a,i,k-1)
       call quickSortP_r8(a,k,j)
    end if
  end subroutine quickSortP_r8
  !------------------------------------------------------!
  subroutine persort_r8(n,a,p)
    integer*4,intent(in)::n
    real*8,intent(inout)::a(n)
    integer*4,intent(out),target::p(n)
    integer*4 i
    do i=1,n
       p(i)=i
    end do
    permutation=>p
    call quickSortP_r8(a,1,n)
    nullify(permutation)
  end subroutine persort_r8
  !------------------------------------------------------!
end module kuty
